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•1 


Summary 


Analytical  and  numerical  methods  were  used  to  investigate  the  flow  induced 
by  a  shock  wave  in  a  shock-tube  channel  containing  air  laden  with  suspended  small 
solid  particles.  Exact  results  are  given  for  the  frozen  and  equilibrium  shock- 
wave  properties  as  a  function  of  diaphragm-pressure  ratio  and  shock-wave  Mach 
numbers.  The  driver  contained  air  at  high  pressure.  A  modified  random-choice 
method  together  with  an  operator- splitting  technique  show  clearly  both  the  decay 
of  a  discontinuous  frozen  shock  wave  and  a  contact  discontinuity  and  the  forma¬ 
tion  of  a  stationary  shock  structure  and  an  effective  contact  front  of  finite 
thickness. 

The  effects  of  particle  diameter,  particle-number  density  and  diaphragm- 
pressure  ratio  on  the  transitional  behaviour  of  the  flow  are  investigated  in 
detail.  The  alteration  of  the  flow  properties  due  to  the  presence  of  particles 
is  discussed  in  detail  and  compared  with  classical  shock-tube  flows. 
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Notation 


equilibrium  speed  of  sound 

frozen  speed  of  sound 

specific  heat  of  solid  material 

gas  specific  heat  at  constant  pressure 

gas  specific  heat  at  constant  volume 

drag  force  acting  on  a  particle 

diameter  of  particles 

gas  thermal  conductivity 

reference  length  [Eq.  (22)] 

mass  of  a  particle 

particle  number  density 

gas  pressure 

rate  of  heat  transfer  to  a  particle 

gas  constant 

gas  temperature 

time  from  •liaphragm  rupture 

gas  velocity 

i  reparation  velocity  of  shock  wave 
velocity  ol'  particles 

space  coordinate  measured  from  diaphragm 

temperature  of  particles 
gas  viscosity 
gas  density 

density  of  solid  particles 
mass  concentration  of  particles 


Dimensionless  Quantities 


CD 

M 

s 

Nu 

P 

P21 

P4l 

Pr 

Re 

T 

'l21 

T31 

u 

U21 

X 

a 

P 

7 

r 

r21 

r3i 

T 

Note: 


=  drag  coefficient  of  a  particle 

=  Valf  °r  Vale 

=  Nusselt  number  of  a  particle 
=  P/Pi 

*  Vpi 
■  *>1 

=  gas  Prandtl  number 
=  Reynolds  number  of  a  particle 
=  l/T,  or  0/T1 

“  VT1 

=  V*1 

=  u/a1  or  v/a^ 

■  vS. 

-  x/i 

=  ratio  of  mass  concentrations  of  particles  to  gas  [Eq.  (19)] 
=  ratio  of  specific  heats  of  two  phases  [Eq.  (l8)] 

=  gas  specific  heat  ratio 

-  specific  heat  ratio  of  an  effective  perfect  gas 
=  p/pj_  or  or/ or^ 

=  P2/Pl 

■  VU/alf) 

Ppi>  Tg^,  U  1  and  can  have  frozen  and  equilibrium  values;  at 
«■  *  -  ibj  1111  r2i  =  P2/pl  = 

vi 


1.  INTRODUCTION 


When  a  gas  carries  a  lot  of  solid  particles,  they  significantly  affect  the 
flow  through  the  transfer  of  momentum  and  heat  from  or  to  the  gas.  Shock  waves 
propagating  in  such  a  dusty  gas  are  characterized  by  a  transition  region  orders 
thicker  than  that  caused  by  viscosity  and  heat  conduction  in  a  pure  gas.  Across 
the  transition  front,  the  interaction  of  the  gas  and  the  particles  leads  to  an 
equilibrium  state  of  the  mixture.  The  structure  of  stationary  shock  waves 
(Refs.  1-5)  and  its  formation  by  a  moving  piston  (Refs.  5,  6)  have  been  studied 
theoretically  on  the  assunption  that  the  formulae  for  the  drag  and  the  rate  of 
heat  transfer  for  a  single  spherical  particle  placed  in  a  steady  flow  can  still 
be  applied  to  the  motion  of  many  particles  contained  in  a  dusty  gas. 

Some  experimental  studies  of  shock  waves  in  a  dusty  gas  inside  a  shock  tube 
were  done  in  order  to  get  some  data  on  the  interaction  of  the  two  phases  (Refs. 
7-10).  Some  of  the  results  showed  an  effective  drag  coefficient,  obtained  from 
the  observation  of  the  acceleration  of  the  particles  behind  the  shock  waves, 
which  differed  appreciably  from  that  for  a  single  particle.  However,  there 
were  many  factors  influencing  the  results  and  a  definitive  conclusion  on  the 
appropriate  drag  coefficient  to  be  used  is  not  available  as  yet. 

Recently,  numerical  analyses  followed  these  experimental  studies.  The 
shock-tube  problem  for  a  dusty  gas  was  solved  numerically  by  Otterman  and 
Levine  (Ref.  ll)  using  the  partiele-in-cell  method.  They  investigated  the 
difference  in  the  transient  flow  field  due  to  different  assunptions  of  the 
drag  and  the  heat  transfer  coefficients.  Outa,  Tajima  and  Morii  (Ref.  10) 
made  a  numerical  analysis  of  the-  penetration  of  a  shock  wave  into  a  dusty-gas 
region  for  conparison  with  their  experimental  results.  Satofuka  and  Tokita 
(Ref.  12)  discussed  the  efficiency  of  various  numerical  techniques  applied  to 
the  case  when  both  sides  of  a  diaphragm  are  filled  with  a  dusty  gas. 

The  work  of  Otterman  and  Levine  provided  a  rough  sketch  of  the  transient 
flow  behaviour  in  a  shock  tribe.  However,  they  treated  cases  of  unusually  high 
mass-loading  ratio  and  their  numerical  results  include  obscure  points  regarding 
the  frozen  shock  front  and  the  delay  of  particle  acceleration. 

At  present,  :  more  extensive  aid  clear  analytical  study  is  required  for 
comparison  with  experimental  results.  In  this  paper,  we  consider  the  classical 
problem  of  the  shock  tube  where  the  driver  contains  high-pressure  air  and  the 
channel  contains  a  dusty  gas,  as  would  be  the  case  in  actual  experiments.  The 
effects  of  the  ratio  of  mass  concentrations,  the  size  of  the  particles  and  the 
diaphragm  pressure  ratio  upon  the  flow  characteristics  are  fully  discussed. 
Working  curves  are  produced  for  the  physical  properties  as  a  function  of  the 
initial  conditions  and  the  diaphragm  pressure  ratio.  The  transitions  through 
the  shock  front  and  contact  region  as  well  as  the  rarefaction  wave  are  studied 
in  detail  as  functions  of  time,  home  of  the  randomness  produced  by  mesh  size 
and  drag  relations  on  physical  properties  are  also  shown. 

The  study  of  shock  waves  in  dusty  gases  provides  a  good  introduction  to 
real-gas  effects.  Frozen  shock  waves,  transition  to  equilibrium  flow,  frozen 
cjid  equilibrium  sound  speeds  and  dispersed  shock  waves  are  all  encountered 
and  well  illustrated  by  analogy  with  real-gas  flows. 


2.  ANALYTICAL  CONSIDERATIONS 


We  consider  the  transient  flow  occurring  after  diaphragm  rupture  in  a 
shock  tube  (Fig.  l) .  In  order  to  formulate  the  motion  of  the  mixture,  we 
must  make  several  assumptions  (Refs.  1-6).  The  gas  is  assured  perfect  and 
its  viscosity  and  heat  conductivity  are  neglected  except  for  the  interaction 
with  the  particles.  The  particles  are  assumed  to  he  spheres  of  uniform  size 
and  their  number  is  so  large  that  the  flow  may  be  treated  as  a  continuum. 

The  volume  occupied  by  the  particles  is  neglected. 

Let  p,  p,  T,  u  he  the  pressure,  density,  temperature  and  velocity  of  the 
gas,  and  u,  9,  v  be  the  mass  concentration,  temperature  and  velocity  of  the 
particles,  respectively.  Using  the  assumptions  stated  above,  we  can  express 
the  equations  of  motion  of  mass,  momentum  and  energy  for  either  the  gas  or 
the  particles  as  follows: 


§e  +h (pu)  ~ 0 


I?  +  (ov)  ~  0 


It  (pu)  +  Sc  ^pu<^  + =  ■  5 D 


3x  m 


a  ^  +  b  (ov2)  =-iD 


St 


!e lp ( cvt  +lu2  u  +  d  ■ /  ~  ~  - 1  2  v  a 


ji  h 


mi  pu(cpt  Jr -i(vD  +  Q) 


It  1  T  (  Cme  +  I  v&  )}  +  h  I  °V  V  V  +  I  y2  )}  =  5  (vD  +  Q) 

The  thermally  perfect  equation  of  state  for  the  gas  is  given  by 


(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


p  =  pRT 


(7) 


The  interaction  of  the  two  phases  is  incorporated  in  the  terms  of  the  drag 
and  the  heat  transfer  on  the  right-hand  sides  of  Eqs.  (3)-(6).  Specification 
of  the  dependence  of  D  and  Q  on  the  flow  quantities  is  needed  to  obtain  a  closed 
set  of  equations.  While  various  formulae  valid  for  a  single  sphere  exist  (Refs. 
3,  4),  attention  should  be  paid  to  the  fact  that  in  the  present  transient  case 
the  Reynolds  number  of  the  particle  takes  on  a  high  value  initially  because  of 
the  large  difference  in  velocity  between  the  two  phases.  We  assume  that 
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D  =  B  d2p(u_v) !u-v |cD 

=  •g  d2p(u-v)  |u-v  |(0.48  +  28  Re-<^*^) 


(8) 


for  the  drag  (Ref.  13)  and 


Q  =  77d(-tC  Pr  'L(T-0)Ku 


«  TTd)iCpPr"1(T-0)(2.O  +  0.6  Pr1/3  Re1//2) 


(9) 


for  the  heat  transfer  (Ref.  l4),  where  Re  is  the  Reynolds  number  based  on  the 
diameter  of  the  particle  and  the  relative  velocity  of  the  particle  to  the  gas 


Re  =  0 |u-v  M fn  (10) 


and  Pr  is  the  Prandtl  number 


Pr  =  |iCp/k  (11) 

The  viscosity  and  the  thermal  conductivity  of  the  gas  are  assumed  constant 
for  stationary  shock  waves  (i.cf.  9),  since  the  change  in  temperature  of  the  gas 
is  small  over  the  relaxation  region.  In  the  present  case,  however,  the  particles 
interact  not  only  with  a  hot  compressed  gas  between  the  shock  front  and  the 
contact  front  but  also  with  a  cold  expanding  gas  behind  the  contact  front.  In 
addition,  the  temperature  of  the  gas  varies  with  time.  Therefore,  the  variations 
of  the  viscosity  and  the  thermal  conductivity  with  temperature  must  be  taken  into 
account.  With  air  (Ref.  15)  inside  the  tube 

,-'.77 

U  =  1.71  x  10  x  \  poise  (12) 

and 

Pr  =  0.75  (13) 


The  high  pressure  air  in  the  driver  obeys  Eqs.  (l),  (3)  and  (5)  with 
a  =  0,  and  Eq.  (7) . 

The  initial  conditions  at  t  =  0  are  as  follows: 


=  £_=!L-=iL=£_  =  i  u  =  v  =  0  for  x  >  0 
1  Pi  T1  ffl  T 


£_  =  £L 
p4  Pit 


u  =  0,  cr  =  0  for  x  <  0 


(14) 
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It  is  assumed  here  that  the  gas  and  the  particles  initially  are  in  equili¬ 
brium.  The  system  of  equations  (l)-(7)  for  seven  unknowns  (p,  p,  T,  u,  c r, 
e,  v)  together  with  the  supplementary  equations  (8) -( 13)  subject  to  the 
initial  conditions  can  be  solved  numerically.  We  apply  the  random-choice 
method  (Refs.  16-19)  to  the  present  problem,  for  it  describes  discontinuities 
clearly  without  the  use  of  implicit  or  explicit  artificial  viscosity.  An 
operator- splitting  technique  makes  its  application  possible  (Refs.  l8,  19). 
The  solution  is  obtained  by  solving  the  two  sets  of  equations  alternately  in 
each  time  step,  that  is,  the  equations  derived  from  Eqs.  (l)-(6)  with  the 
right-hand  sides  omitted  and  the  equations  obtained  by  omitting  the  deriva¬ 
tives  with  respect  to  x. 

The  element  of  the  random-choice  method,  which  solves  the  former  set, 
is  to  determine  by  means  of  random  sampling  the  solution  at  the  middle  point 
between  the  two  adjacent  points  where  the  solution  is  known  at  a  previous 
time.  The  initial  condition  in  this  calculation  is  taken  as  a  step-like 
discontinuity.  The  procedure  for  the  gas  phase  has  been  written  in  a  few 
papers  (Refs.  17-19)  and  is  omitted  in  this  paper. 

On  the  other  hand,  a  difficulty  arises  in  the  treatment  of  the  equations 
for  the  particles  in  this  stage,  which  are  reduced  to 


ckr  b 
3t  +  35 


(ov) 
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8v 

St 


+  V 


dv 

35 
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at 

3t 


+  V 


So 


=  0 


(15) 


A  discontinuous  initial  condition  would  result  in  a  multivalued  solution  for 
the  velocity  v.  The  flow  variables  of  the  particles  are  continuous  and  we 
assume  the  initial  condition  of  linear  distribution  only  for  the  particles  in 
order  to  avoid  this  difficulty.  Consider  two  adjacent  points  with  distance 
Ax  apart,  where  the  solution  is  known.  We  take  the  following  initial  condi¬ 
tions: 


9  "  (9r  "  V  h  +  \  (0r  +  V 

a  =  (°r  -  °j)  h  +i  (0r  +  ^ 


(16) 


where  the  subscripts  r  and  l  refer  to  the  two  points  x  =  ±  §■  Ax  (for  convenience), 
respectively.  The  solution  of  Eq.  (15)  subject  to  the  condition  (l6)  is  given  by 
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Random  sampling  of  the  solution  between  the  two  points  determines  the  values 
of  flow  variables  at  the  middle  point  a  little  time  later.  Two  repetitions 
of  this  procedure  bring  about  a  time  evolution  by  a  time  step,  say,  At.  The 
time  step  must  be  appropriately  small  so  that  the  so-called  Courant-Friedrichs- 
Lewy  condition  is  satisfied,  otherwise  the  result  loses  physical  meaning.  Al¬ 
though  most  points  are  treated  in  this  way,  the  point  bordering  the  region  of 
the  dusuy  gas  is  an  exception.  Since  the  diffusion  of  concentration  is  not 

oonsj  •’ . ■'  m- -t,  *-  a  l4  s'r-.ti  : :  initial  condition  at  this  point  so 

as  co  got  a  clear  oouuuary  of  the  particle  cloud. 

The  other  portion  of  the  operator- splitting  technique  is  accomplished 
by  integration  of  the  equations  with  respect  to  time.  This  procedure  must  be 
carefully  done  for  moving  discontinuities,  because  it  may  produce  a  significant 
error  of  the  first  order.  Consider  the  case  when  a  discontinuous  surface 
passes  across  a  point  in  a  tine  of  At.  The  surface  should  be  regarded  to  have 
moved  over  the  distance  Ax  between  the  middle  points  on  both  sides  of  the  point. 
Therefore,  the  contribution  to  the  integration  in  At  at  the  point  must  be  con- 
sidered  separately  for  the  two  halves  of  |  At. 


3.  RECULTS  AND  jjlXbSait 

An  assignment  of  several  parameters  is  needed  for  the  numerical  calculations. 
We  take  air  (7  =  l.U)  as  the  gas  in  both  the  driver  and  the  driven  sections.  The 
is  assumed  to  o«  at  room  temperature  (Tq  =  Ti+ )  before  diaphragm  rupture.  The 
initial  pressure  pj_  in  the  driven  section  is  taken  to  be  one  atmosphere.  Further, 
we  restrict  ourselves  to  the  case  where  the  ratio  of  specific  heats  of  the  two 
phases 

P  =  C  /C  (18) 

m'  v 

is  unity.  The  regaining  factors,  on  which  flow  behaviours  depend,  are  the  ratio 
of  mass  concentration 

a  =  (19) 

the  diameter  of  the  particles  d  and  the  diaphragm  pressure  ratio 

*41  =  P4/P1 
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The  results  for  the  case  of  a  =  1,  d  =  10  Jim  and  P4i  =  10  are  presented 
in  Figs.  2-6.  Normalized  quantities  are  used  in  all  the  figures.  Flow  quan¬ 
tities,  except  the  velocity,  are  referred  to  corresponding  values  in  the 
stationary  region  (1)  ahead  of  the  shock  wave.  The  velocity  is  normalized 
by  the  speed  of  sound  in  the  gas  phase,  i.e.,  the  so-called  frozen  speed  of 
sound: 


alf  =  'Iw1/p1 

The  distance  x  is  also  measured  in  units  of 


(21) 


8m 

-2 

TTP-jd 


k 

3 


(22) 


and  the  time  t  in  i/aif.  [Note  the  value  of  JL  can  be  obtained  from  Eq.  (1)]. 
Numerical  calculations  have  been  done  in  the  manner  stated  in  the  previous 
section.  A  mesh  size  Ax  of  0.1  is  used  together  with  time  step  At  as  large  as 
the  Courant-Friedrichs-Lewy  condition  may  allow. 


Flow  structures  after  a  small  time  has  elapsed  after  the  diaphragm  rupture 
are  shown  in  Fig.  2,  for  the  case  of  P4l  =  10,  a  =  1,  d  =  10  pm  for  a  time  t  -  h . 
Under  these  conditions  with  pq  =  1  atm  and  Tq  =  300  K,  using  glass  beads,  the 
number  density  to  provide  a  =  1  would  be  0.9^  x  lO^/cc;  l  =  2.72  cm;  aqf  =  350 
m/s  and  t  =  4  s  3.11  x  10“^  sec  (see  Table  1  for  further  details).  The  dashed 
lines  are  for  the  classical  shock-tube  problem.  The  solid  lines  show  the  frozen 
shock  wave  as  an  abrupt  change  followed  by  a  gradual  approach  to  equilibrium. 

In  this  case  the  flow  is  far  from  equilibrium.  Note  in  Fig.  2(b)  the  particle 
concentration  rises  gradually  though  the  shock  front,  reaches  a  maximum  at  the 
contact  front  and  then  drops  to  zero  at  the  driver  gas.  It  should  be  noted  that 
the  rarefaction  is  weaker  than  for  the  gas  case  only.  In  addition  there  are 
spurious  numerical  oscillations  near  the  tail  of  the  wave  in  all  flow  properties. 
Only  the  gas  responds  to  the  abrupt  change  at  the  instant  of  diaphragm  rupture, 
while  the  particles  cannot  folow  any  sudden  change.  After  the  frozen  shock 
front  has  passed  by,  the  velocity  and  tenperature  of  the  particles  are  raised 
gradually  through  interaction  with  the  gas.  On  the  other  hand,  the  gas  is 
decelerated  and  loses  energy.  A  comparison  of  the  results  with  the  solution 
for  a  pure  gas  exhibits  a  decay  of  the  frozen  shock  front  due  to  this  inter¬ 
action.  The  velocity  of  propagation  of  the  shock  wave  also  diminishes  as 
can  be  seen  by  looking  at  the  values  of  X  for  both  cases.  The  deceleration 
of  the  gas  results  in  a  compression  and  its  pressure  away  from  the  frozen 
chock  front  attains  a  higher  value  than  in  the  case  of  a  pure  gas.  The 
rarefaction  wave  weakens  as  a  result.  Some  particles  cross  the  frozen  contact 
surface  and  drop  to  zero  in  the  cold  gas.  Thus  the  temperature  of  these  particles 
drops.  However,  the  gas  temperature  [see  Fig.  2(c)]  also  drops  but  not  as  low 
as  for  the  frozen  flow  since  the  particles  act  as  a  heat  reservoir. 

Subsequent  transitional  behaviours  of  the  physical  quantities  are  shown  in 
Figs.  3-5  for  increasing  time  t.  Both  the  velocities  and  temperatures  of  the  gas 
and  the  particles  behind  the  frozen  shock  front  approach  each  other  with  time  and 
a  new  uniform  region  in  equilibrium  forms  some  distance  behind  the  discontinuous 
frozen  shock  front.  An  almost  stationary  shock  structure  of  finite  thickness 
can  be  seen  in  Fig.  5.  For  our  case  of  Pp  =  2.5  g/cm3,  this  state  is  attained 
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after  2.49  x  10"3  sec  and  the  thickness  of  the  shock  transition  is  about  34  cm. 
If,  for  example,  Pp  is  made  smaller  for  a  =  1,  the  particle  number  density  will 
increase  and  the  transition  length  will  decrease. 

It  is  very  useful  to  assume  that  everywhere  the  particles  reach  the  equili¬ 
brium^  flow  limit  of  the  gas  velocity  and  temperature.  This  must  occur  because^ 
the  flow  except  for  relaxation  regions  of  finite  length  must  approach  this  limit 
after  a  sufficient  time  elapses.  The  mixture  in  this  limit  behaves  effectively 
as  a  perfect  gas  with  the  specific  heat  ratio  given  by  (Ref.  4) 

y^  =  ( y  +  Qp)/(l  +  Cd£)  (23) 

and  the  speed  of  sound  by 

(24) 


The  shock-tube  solution  for  this  equilibrium-flow  limit  is  also  shown  in  Fig.  5 
10 r  comparison.  The  agreement  of  the  equilibrium  flow  limit  and  the  numerical 
results  is  very  good  indeed.  Besides  the  relaxation  region  of  the  shock  wave, 
the  contact  surface  also  has  a  structure  of  finite  thickness.  The  particles  in 
this  layer  are  in  equilibrium  with  the  gas.  The  structure  of  the  effective 
contact  surface  must  depend  on  how  the  particles  have  interacted  with  the  gas. 

A  wave  diagram  is  shown  in  Fig.  6.  The  paths  of  the  discontinuous  frozen 
and  equilibrium  shock  fronts  are  shown.  For  comparison  the  present  numerical 
results  are  also  shown.  As  noted  on  Fig.  6,  the  present  data  shows  the  relaxa¬ 
tion  distance  between  the  present  frozen  shock  and  the  idealized  equilibrium 
gas-particle  front.  It  shows  that  the  present  frozen  shock  moves  at  a  velocity 
closer  to  the  idealized  equilibrium  shock  front  than  the  perfect  gas  frozen 
shock  wave  after  a  sufficient  time  has  elapsed.  At  times  close  to  zero  it 
would  move  closer  to  the  perfect  gas  shock. 

The  path  of  the  contact  surface  has  been  made  possible  by  choosing  a 
point  on  its  transition  region  such  that  the  temperature  equals  the  original 
gas  temperature  Tq .  Consequently,  the  results  by  definition  must  lie  on  our 
equilibrium  contact  front  path  for  convenience  of  illustration.  The  results 
for  the  tail  of  the  rarefaction  wave  were  obtained  by  extrapolating  the 
straight  lines  of  the  pressure  curves  in  regions  (2)  and  (3)»  the  uniform 
states,  until  it  hit  the  rarefaction  wave  profile.  It  can  be  seen  that  the 
agreement  with  the  idealized  equilibrium  value  is  very  good  by  large.  Note 
that  the  tail  of  this  wave  is  weaker  than  the  perfect  gas  tail.  Of  course 
the  head  of  the  rarefaction  wave  is  identical  for  all  three  flows  as  shown. 

The  changes  in  temperature  and  pressure  of  the  gas  immediately  behind 
the  frozen  shock  front  as  it  moves  are  shown  in  Fig.  7*  They  start  from  values 
for  a  pure  gas  and  approach  finally  the  values  calculated  from  the  shock  speed 
attained  in  the  idealized  equilibrium- flow  limit. 
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Small  random  disturbances  in  the  numerical  results  are  characteristic 
of  the  use  of  the  random-choice  method.  Taking  a  smaller  mesh  size  makes 
the  disturbances  smaller.  (We  confirmed  this  by  comparing  the  results  with 
those  for  a  half  and  quarter- size  mesh  smaller,  although  only  for  a  short 
time  owing  to  the  increased  cost  of  the  computations  -  see  Figs.  8  and  9>) 

It  should  be  noted  that  the  thermal  equilibration  between  the  two  phases 
is  achieved  faster  than  the  equilibration  of  velocity  (see  Figs.  2c,  d  and 
3c,  d) .  This  fact  is  also  reflected  in  that  longer  length  from  the  discontinuous 
shock  front  is  needed  for  velocity  equilibration  than  for  temperature  equilibra¬ 
tion  (see  Figs.  5c  and  51). 

Next,  we  investigate  the  influence  of  the  size  of  the  particles  on  the 
transition  of  the  How.  The  results  for  the  cases  of  d  =  20  pm  (see  Figs. 

10-12)  and  d  =  40  pm  (see  Figs.  13-15)  with  other  parameters  unchanged,  are 
shown  in  Figs.  10-15,  respectively.  It  is  interesting  to  find  almost  similar 
flow  fields  at  t  =  )  in  the  three  cases  of  10-40  pm.  Note  that  the  actual 
times  and  distances  aie  different  for  different  particle  diameters,  being 
proportional  to  the  cJ.ar.eter  according  to  Eq.  (22).  This  is  nearly  true  for 
later  times  as  well .  This  similar:! nay  be  minted  to  the  constancy  of  the 
drag  coefficient  of  the  particles  feu  large  Reynolds  numbers.  The  rate  of 

heat  transfer  is  much  greater  for  smaller  particles.  For  example,  a  comparison 

of  Figs.  2c  and  13c  soon  show  that  the  10  pm  particle  differs  greatly  from  the 
frozen  value  but  the  to  pm  particle  is  very  close  to  the  frozen  value. 

The  results  at  t  =  32  may  be  compared  with  the  structure  of  a  stationary 

shock  wave.  It  can  be  seen  that  the  stationary  shock  values  are  almost' 

achieved.  The  length  of  the  relaxation  region  is  naturally  longer  for  larger 
particles  on  account  of  their  larger  inertia.  Alternatively,  the  relaxation 
length  depends  on  the  diameter  d  [Eq.  (22)].  However  it  is  found  to  lie  between 
d  and  d  .  As  seen  from  Figs.  7,  12  and  15,  the  time  required  for  the  stationary 
shock  wave  to  form  also  depends  on  the  diameter  of  the  particles.  However,  it 
is  not  a  linear  relation.  The  time  t  for  the  20  pm  particle  to  achieve  equili¬ 
brium  is  much  less  them  a  factor  of  two  for  the  10  pm  particles.  Nevertheless 
the  actual  time  is  longer  than  two  but  less  than  four. 

Calculi.. fc; o.i;  j,vo  been  done  for  cases  of  different  ratios  of  mass  concen¬ 
trations  with  the  remaining  parameters  the  same  as  the  first  case.  The  results 
for  a  =  0.1  and  a  =  2.0  are  shown  in  Figs.  l6-l8  and  19-21  respectively.  It 
can  be  seen  from  comparison  of  the  flow  fields  at  t  =  4  that  a  larger  mass 
concentration  of  particles  causes  in  the  gas  phase  larger  deviations  from  the 
frozen  flow.  However,  the  final  quasi-equilibrium  state  is  accomplished  more 
quickly  when  a  is  small.  For  example,  compare  median  curves  through  Figs.  18, 

7  and  21  and  it  will  be  clear  that  equilibrium  is  faster  for  the  small  loading 
ratio  a  =  0.1.  The  effect  of  the  mass  concentration  of  particles  on  the  thick¬ 
ness  of  the  stationary  shock  wave  is  not  so  apparent.  In  the  case  when  Q  =  2.0, 
the  propagation  of  the  shock  wave  becomes  so  slow  that  the  discontinuous  Jump 
in  the  gas  phase  cannot  be  supported  (see  Fig.  20).  The  transition  to  the 
fully  dispersed  shock  wave  is  characteristic  of  this  case.  The  contact  surface 
consists  of  a  region  of  dusty  gas  of  finite  thickness  which  is  followed  by  a 
discontinuity  in  the  gas  phase  (see  Figs.  20b  and  20c). 

The  difference  in  transitional  behaviour  due  to  the  strength  of  the  diaphragm 
pressure  ratio  is  studied  next.  The  results  for  the  cases  of  ptp  =  5 .0  and  pj4p=  20 
are  shown  in  figs.  22-24  and  25-27,  respectively.  The  length  of  the  relaxation 
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region  of  the  stationary  shock  wave  forming  after  a  sufficient  time  is  larger 
for  weaker  shock  waves  (see  Figs.  23  and  26),  i.e.,  smaller  diaphragm  pressure 
ratio.  When  the  diaphragm  pressure  ratio  is  lower  than  a  critical  value,  the 
shock  wave  is  weak  and  the  decay  due  to  absorption  of  energy  by  the  particles 
is  so  large  that  the  shock  wave  becomes  fully  dispersed  ,  that  is,  the  frozen 
shock  front  disappears.  In  fact,  the  case  of  pi^  =5.0  lies  in  this  range. 
Comparison  of  the  changes  in  the  pressure  jump  at  the  frozen  shock  front  with 
time  clarifies  that  the  stationary  shock  wave  forms  also  in  a  shorter  time  for 
a  higher  diaphragm  pressure  ratio  (see  Figs.  7,  24  and  27). 

Finally,  several  figures  (Figs.  28-34)  are  presented  illustrating  the 
effect  of  the  existence  of  particles  upon  the  uniform  states  between  the  shock 
wave,  the  contact  surface  and  the  rarefaction  wave  for  the  idealized  gas-particle 
equilibrium  flow  limit.  Flow  quantities  are  given  by  the  exact  classical  shock- 
tube  solution  (Ref.  20)  (see  Appendix  B)  for  the  effective  perfect  gas  based 
on  yc  and  alf  (see  Table  2).  Variations  of  shock  Mach  numbers  based  on  the 
frozen  or  the  equilibrium  speed  of  sound  with  the  diaphragm  pressure  ratio  axe 
shown  in  Fig.  23  for  values  of  cc  over  the  range  0  <  a  <  2.  The  particles  reduce 
the  velocity  of  propagation  of  a  shock  wave  and  since  the  speed  of  sound  of  the 
gas  is  fixed  the  Mach  number  Ms  falls  with  increasing  a  for  a  fixed  pltp.  The 
reason  lies  in  the  fact  that  the  particles  absorb  kinetic  and  thermal  energy, 
however,  the  equilibrium  speed  of  sound  becomes  much  smaller  and  as  a  result 
the  effective  shock  Mach  number  increases  with  a  for  a  fixed  PM.  The  shock 
speed  is  even  less  than  the  frozen  speed  of  sound  if  the  amount  of  particles 
is  sufficiently  large.  Then,  the  shock  wave  becomes  dispersed  without  a 
discontinuous  frozen  front.  The  region  below  the  dashed  line  Ms  =  1,  is  that 
of  dispersed  shock  waves  over  the  pi^  range.  The  variations  of  flow  quantities 
behind  the  shock  wave  with  diaphragm  pressure  ratio  are  illustrated  in  Figs. 
29-32.  The  flow  quantities  of  the  gas  immediately  behind  the  frozen  shock  front 
are  also  plotted.  The  temperature  and  the  velocity  reduce  as  cc  takes  on  larger 
value.  On  the  contrary,  the  equilibrium  pressure  and  density  increase  owing  to 
compression.  It  should  be  noted  that  at  equilibrium  the  frozen  values  of  the 
pressure  and  the  density  of  the  gas  decreases  with  particle  concentration.  The 
density  and  the  temperature  between  the  contact  surface  and  the  rarefaction 
wave  are  plotted  in  Figs.  33  and  34.  Tbe  presence  of  particles  bring  about 
compressive  effects  on  the  flow  in  this  region,  such  that  these  values  increase 
with  particle  concentration  for  a  given  diaphragm  pressure  ratio. 

4.  COIJCLUSIOHo 

Flow  properties  in  a  shock  tube,  in  which  many  solid  particles  are  sus¬ 
pended  in  the  driven  section,  were  analyzed  numerically.  Use  was  made  of  the 
r an lom- choice  method  modified  so  as  to  be  applicable  to  a  dusty  gas  together 
with  an  operator  splitting  technique. 

The  particles  remove  momentum  and  energy  from  the  gas  behind  the  shock 
wave.  As  a  result,  the  strength  of  the  discontinuous  frozen  shock  wave  in 
the  gas  phase  decays.  If  the  number  of  particles  is  sufficiently  large  or 
the  shock  is  fairly  weak,  the  frozen  shock  decays  to  a  Mach  wave  and  the 
shock  wave  becomes  fully  dispersed.  The  deceleration  of  the  gas,  while  the 
particles  accelerate  behind  the  shock  wave,  produces  a  compression  and  the 
equilibrium  pressure  becomes  higher  than  for  the  case  of  a  pure  gas  at  some 
distance  from  the  front.  In  this  manner  a  thick  stationary  struct\xe  arises, 
where  the  particles  finally  reach  the  same  velocity  and  temperature  as  the 
gas.  On  the  other  side  of  the  contact  front  the  compressive  effect  weakens 
the  strength  of  the  rarefaction  wave. 


Some  particles  remain  in  the  contact  region,  such  that  they  interact 
with  the  hot  gas  in  front  of  it  and  with  the  cold  gas  behind  it.  Thus,  the 
sharp  discontinuities  in  temperature  and  density  of  the  gas  become  transi¬ 
tional  and  a  contact  region  of  finite  thickness  appears.  For  the  cases  con¬ 
sidered  here,  with  larger  drag  values,  it  is  possible  for  the  particles  to 
vanish  abruptly  at  the  cold  side  of  the  contact  region.  As  a  consequence, 
there  is  a  sudden  discontinuity  in  temperature  and  density  typical  of  a  contact 
surface  (see  Fig.  20). 

The  influence  of  particle  diameter,  particle-number  density  and  diaphragm 
pressure  upon  the  transient  properties  of  the  flow  have  been  studied  with  the 
following  conclusions .  Particles  of  large  size  increase  the  time  required  for 
the  flow  to  be  in  quasi-equilibrium  and  also  increase  relaxation  length  or  time 
to  the  final  stationary  shock  wave.  The  degree  of  increase  varies  with  particle 
diameter  d  and  lies  between  d  and  d^ .  When  the  particle  number  density  is 
large,  variations  in  the  flow  quantities  occur  quickly  but  the  final  equilibrium 
flow  values  are  reached  after  a  longer  time.  Strong  shock  waves  have  much 
shorter  relaxation  lengths. 

The  flow  quantities  in  the  equilibrium-flow  limit  were  calculated  from 
exact  shock-tube  relations.  The  speed  of  sound  and  the  specific-heat  ratio 
of  this  effective  perfect  gas  take  on  fairly  low  values  owing  to  the  increase 
of  the  effective  molecular  weight.  The  increase  in  pressure  and  density  of 
the  dust-laden  gas  from  its  frozen  values  at  the  shock  wave  can  be  quite  large. 
The  temperature  and  velocity  on  the  other  hand  decrease. 

The  results  obtained  in  this  paper  have  been  based,  in  particular,  on  the 
assumption  that  the  drag  force  and  the  rate  of  heat  transfer  to  the  particles 
are  given  by  Eqs.  (8)  and  (9)  (see  Appendix  C  for  other  assumptions) .  however, 
there  are  many  causes  in  practice  to  make  these  assumptions  questionable.  Among 
these  are  non-spherical  shape  of  particles,  variation  in  local  distribution  of 
particle- size,  interaction  between  particles,  rotation  of  particles  and  electro¬ 
static  forces.  Nevertheless,  the  quantitative  and  qualitative  nature  of  the 
flow  in  a  dusty-gas  shock  tube  has  been  made  clear  in  our  study. 

Undoubtedly  improvements  will  be  made  in  the  future  to  take  into  account 
some  of  the  above  non-ideal  properties  of  the  particles.  It  is  possible  that 
even  volume  and  partial -pressure  effects  will  be  considered  in  future  for 
comparison  with  the  present  study. 
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APPENDIX  A 


LIST  OF  NUMERICAL  PROGRAMS 


In  each  time  step,  four  calculations  are  done  in  order:  integration  of 
inhomogeneous  equations  for  §  At,  two  applications  of  random-choice  procedure 
to  homogeneous  equations  and  integration  of  inhomogeneous  equations  for  \  At . 
Godunov's  iterative  procedure  (Refs.  17-19)  is  used  in  solving  the  Riemann 
problem  for  the  gas. 
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APPENDIX  B 


FROZEN  AND  EQUILIBRIUM  FLOWS 


The  frozen-flow  values  at  the  instant  of  diaphragm  rupture  when  the 
particles  have  no  effect  can  readily  be  found  from  the  relations  given  in 
Ref.  20.  If  it  is  now  assumed  that  the  velocity  and  the  tenperatuxe  of  the 
particles  are  the  same  as  those  of  the  gas  everywhere  and  we  neglect  the 
transition  thicknesses  of  the  shock  wave  and  the  contact  region,  the  equi¬ 
librium-flow  limit  is  readily  found.  In  this  limit,  the  pressure  in  the 
uniform  region  behind  the  shock  wave  p2  must  satisfy  the  shock-tube  equation 
(Ref.  20): 

(74  -  l)(a1/a4)(P2/P1  “  !) 

^2 y1  +  {y1  +  1)(p2/p1  -  l) 


Pi '  px  1. 


_2JU_ 

|  V1  (Bl) 


■n.u  ■  e  uua  ,l.:  -  r-  now  given  bj  ana  (24)  as  7^  and  ap,  respectively. 

Once  P22  is  known,  the  other  flow  quantities  are  obtained  from  the  Rankine- 
Hugoniot  relations  as  follows  (see  Ref.  20): 


p2  _  1  +  (7i  +  1V(r1  -  1)  x  (Pg/p-^ 
(7X  +  1)/(71  -  1)  +  P2/Pq 


(B2) 


u 


wiiej’e  u.,  is  the  velocity  oi 

kJ  0 


-  1)/(71  -  1)  +  P2/P1 

+  1)/(7i  -  1)  X  (P2/Pi) 


1  1  -  1  7p  +  1  P2 
'/  ?-71  *  P^ 

•opagation  of  the  shock  wave. 
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Jhe  tempera:, are  and  tiie  density  in  the  uniform  region  between  the  contact 
surface  and  the  rarefaction  wave  are  given  by  the  isentropic  relations  (Ref.  20) 
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1  *7*1 
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P3  _ /P3  y4 
P4  \  P4  ) 

Once  us  is  obtained,  we  can  calculate  the  properties  of  the  discontinuous 
frozen  shock  wave  at  the  head  of  the  shock  front.  If  u3  is  smaller  than  the 
frozen  speed  of  sound  aIf,  the  frozen  shock  wave  reduces  to  a  Mach  wave  and 
it  is  fully  dispersed. 

For  us  larger  than  alf ,  the  pressure  immediately  behind  the  frozen  shock 
front  is  obtained  from  Eq.  (b5)  by  taking  aif  and  7  as  a^  and  7^  respectively. 
Other  quantities  at  this  position  are  calculated  from  the  foregoing  Rankine- 
Hugoniot  relations  for  a  pure  gas. 

To  summarize,  in  Fig.  26a  for  example,  the  equilibrium  value  of  pressure 
can  be  calculated  at  once  from  the  shock- tube  relations  for  the  given  initial 
conditions.  Knowing  the  equilibrium  Mach  number,  the  frozen  part  of  the  shock 
front  can  also  be  immediately  calculated  from  the  shock- tube  equations  (since 
the  entire  equilibrium  phenomenon  and  the  frozen  shock  all  move  at  the  same 
Macv  "’irber' .  Tic  transition  from  the  frozen  shock  wave  to  equilibrium, 
however,  must  be  obtained  numerically.  On  this  basis  all  the  curves  from 
Fig.  28  to  32  were  constructed  to  provide  exact  solutions  useful  for  the 
experimenter . 
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APPENDIX  C 


NONEQUILIBRIUM  SHOCK-TUBE  PROFILES  USING  STOKES 1  DRAG  LAW 


Figures  35-37  are  included  to  illustrate  how  the  Stokes'  drag  law  Cp  =  Re/24 
together  with  Nu  =  2  can  give  erroneous  results.  For  example,  the  shock  transi¬ 
tion  profiles  are  unduly  long  owing  to  the  very  small  drag  coefficient  (compare 

Fig.  36a  with  Fig.  5a)-  The  same  is  true  of  the  contact  front  transition  (compare 

Fig.  36c  with  Fig.  5c).  The  change  in  transition  is  slow  and  the  oscillations 

are  very  much  reduced  ( compare  Fig .  7  with  Fig .  37) . 
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Table  1 


Some  Reference  Lengths,  Times  and  Number  Penalties 

o  ^  Pp  ,  ,  i 

i  =  *  —  d,  t  = -  t 

3  P-^  alf 

o  o  o 

p  =  2.5  g/enr  (typical  of  crown  glass),  P^  =  1.23x10  J  g/cur ,  aif  =  350  m/sec 


d  (nm) 

10 

20 

40 

i  (cm) 

•  2.72 

5.44 

10.9 

i/alf  (sec) 

0.78  x  10“ 4 

1.56  x  10-4 

-4 

3-11  x  10 

t  (sec)  for  t  =  4 

-4 

3.11  x  10 

-k 

6.23  x  10 

1.25  x  10"3 

t  (sec)  for  t  =  32 

2.49  x  10-3 

4.98  x  10~3 

9.96  x  10”3 

n  (cm-3)  for  a  =  1.0 

0.94  x  106 

1.17  x  105 

1.46  x  104 

Table  2 

Some  Properties  of  Idealized  Equilibrium  Gas-Particle  Mixture 

7  +  Off  /  7  +  Oft ^1 

7e  1+00’  aie  -7  (1  +  a)(l  +  ap)  '  p1 

0=1  (typical  of  glass),  T.^  =  300°K 
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FIG.  2  FLOW  QUANTITIES  AT  t  =  4  (a  = 
-  GAS.  -  PARTICLES, 


(b)  Mass  Concentration 


(d)  Velocity 
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1>  =  10,  d  =  10  jim) 

. -  FROZEN  FLOW 
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FIG.  3  FLOW  QUANTITIES  AP  t  -  8  (o  ■  1,  -  10,  d  -  10  jm) 

-  GAS,  -  PARTICLES. 
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FIG.  5  FLOW  QUANTITIES  AT  t  =  32  (<*  =  1,  =  10,  d  =  10  jm) 

-  GAS,  -  PARTICLES, . EQUILIBRIUM  FLOW 
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FIG.  7 
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FIG .  8  VARIATION  WITH  TIME  OF  TEMPERATURE  AND  PRESSURE  OF  GAS 
JUST  BEHIND  DISCONTINUOUS  FROZEN  SHOCK  FRONT  (a  =  1, 
P^  =  10,  d  =  10  um,  Ax  =  0.05) 

.  EXPECTED  FINAL  VALUE. 


Fit; .  10  FLOW  QUANTITIES  AT  t  =  4  (a  =  1,  ?hl  ~  1°>  d  =  20  |am) 
-  GAS.  -  PARTICLES,  - FROZEN  FLOW. 


o.co 


FIG.  12  VARIATIONS  WITH  TIME  OF  TEMPERATURE  AND  PRESSURE  OF  GAS 
JUST  BEHIND  DISCOOTINUOUS  FROZEN  SHOCK  FRONT  (a  =  1, 

=  10,  d  =  20  jam) 

.  EXPECTED  FINAL  VALUE. 


00 '0 


PARTICLES , 


EQUILIBRIUM  FLOW 


lu.  17  FLOW  QUANTITIES  AT  t  =  32  («  =  0.4,  P^  =  10,  d  =  10  | 
-  GAS,  -  PARTICLES, . EQUILIBRIUM  FLOW 


17  (continued)  flow  QUANTITIES  AT  t  =  32  (a  =  O.U,  P)4l 

d  =  10  p.m)  . 

-  GAS,  -  PARTICLES, . EQUILIBRIUM  FLOW 


FIG  .  18  VARIATIONS  WITH  TIME  OF  TEMPERATURE  AND  PRESSURE  OF  G 
JUST  BEHIND  DISCONTINUOUS  FROZEN  SHOCK  FRONT  (a  =  0.^ 
PI4 1  =  10 ,  d  =  10  nm)  . 

.  EXPECTED  FINAL  VALUE. 
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FIG.  19  FLOW  QUANTITIES  AT  t  =  4  (a  =  2 ,  Pl^  =  10,  d  =  10  urn) 
-  GAS,  -  PARTICLES, - FROEEN  FLOW. 


FIG.  20  (CONTINUED)  FLOW  QUANTITIES  AT 
d  =  10  Jim)  . 

-  GAS,  - 


PARTICLES , 


0.  '40 


FLOW  QUANTITIES  AT  t  =  32  (a  =  l) ,  =  20,  d  =  10  (im) 

-  GAS. - PARTICLES.  . EQUILIBRIUM  FLOW. 


FIG.  26  (CONTINUED)  FLOW  QUANTITIES  AT  T  =  32  (a  =  1,  =  20 

d  =  10  nm)  • 

-  GAS,  - 


PARTICLES , 


EQUILIBRIUM  FLOW. 


FIG.  28  VARIATION  OF  SHOCK  WAVE  MACH  NUMBER  Ms  WITH  DIAPHRAGM  PRESSURE  RATIO  P^. 

-  Ms  9ASED  ON  FROZEN  SPEED  OF  SOUND  a1  ,  -  Ms  BASED  ON  EQUILIBRIUM 

SPEED  OF  SOUND  ai  .  1 


VARIATION  OF  SHOCK  PRESSURE  RATIO  P2i  WITH  DIAPHRAGM  PRESSURE  RATIO  PL,].. 
-  FROEEN  VALUE,  -  EQUILIBRIUM  VALUE. 


VARIATION  OF  DENSITY  RATIO  r21  WITH  DIAPHRAGM  PRESSURE  RATIO  Pj4l 
-  FROZEN  VALUE,  -  EQUILIBRIUM  VALUE. 


:0N  OF  TEMPERATURE  RATIO  T2i  WITH  DIAPHRAGM  PRESSURE  RATIO  P)4l. 
FROZEN  VALUE,  -  EQUILIBRIUM  VALUE. 


VARIATION  OF  VELOCITY  U2  WITH  DIAPHRAGM  PRESSURE  RATIO  P)(1. 
-  FROZEN  VALUE,  -  EQUILIBRIUM  VALUE. 


02  01 


VARIATION  01  DENSITY  RATIO  Ton  WITH  DIAPHRAGM  PRESSURE  RATIO  Pj 


cO  1.02.0 


FIG.  34  VARIATION  OF  TEMPERATURE  RATIO  T3]_  WITH  DIAPHRAGM  PRESSURE  RATIO  P ^ 
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LOW  OLA -TITIEO  AT  =  32  (Cd  =  2U/Re,  Nu 
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FIG.  37  VARIATIONS  WITH  TINE  OF  TEMPERATURE  AND  PRESSURE  OF  GAS 
JUST  BEHIND  DISCONTINUOUS  FROZEN  SHOCK  FRONT  (CD  =  2VRe 
Nu  =  2,  a  =  1,  PUi  =  LO,  d  =  10  nm) . 

.  EXPECTED  FINAL  VALUE. 


